{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 18\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from the previous chapter\n",
    "\n",
    "Read the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('data/glucose_insulin.csv', index_col='time');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolate the insulin data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = interpolate(data.insulin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The glucose minimal model\n",
    "\n",
    "I'll cheat by starting with parameters that fit the data roughly; then we'll see how to improve them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = Params(G0 = 290,\n",
    "                k1 = 0.03,\n",
    "                k2 = 0.02,\n",
    "                k3 = 1e-05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a version of `make_system` that takes the parameters and data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params, data):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "    \n",
    "    params: sequence of G0, k1, k2, k3\n",
    "    data: DataFrame with `glucose` and `insulin`\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    G0, k1, k2, k3 = params\n",
    "    \n",
    "    Gb = data.glucose[0]\n",
    "    Ib = data.insulin[0]\n",
    "    I = interpolate(data.insulin)\n",
    "    \n",
    "    t_0 = get_first_label(data)\n",
    "    t_end = get_last_label(data)\n",
    "\n",
    "    init = State(G=G0, X=0)\n",
    "    \n",
    "    return System(params,\n",
    "                  init=init, Gb=Gb, Ib=Ib, I=I,\n",
    "                  t_0=t_0, t_end=t_end, dt=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the update function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Updates the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: State object\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    k1, k2, k3 = system.k1, system.k2, system.k3 \n",
    "    I, Ib, Gb = system.I, system.Ib, system.Gb\n",
    "    dt = system.dt\n",
    "        \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    G += dGdt * dt\n",
    "    X += dXdt * dt\n",
    "\n",
    "    return State(G=G, X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running the simulation, it is always a good idea to test the update function using the initial conditions.  In this case we can veryify that the results are at least qualitatively correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "update_func(system.init, system.t_0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now `run_simulation` is pretty much the same as it always is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "        \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    init = system.init\n",
    "    t_0, t_end, dt = system.t_0, system.t_end, system.dt\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how we run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_simulation(system, update_func);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are in a `TimeFrame` object with one column per state variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following plot shows the results of the simulation along with the actual glucose data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot(2, 1, 1)\n",
    "\n",
    "plot(results.G, 'b-', label='simulation')\n",
    "plot(data.glucose, 'bo', label='glucose data')\n",
    "decorate(ylabel='Concentration (mg/dL)')\n",
    "\n",
    "subplot(2, 1, 2)\n",
    "\n",
    "plot(results.X, 'C1', label='remote insulin')\n",
    "\n",
    "decorate(xlabel='Time (min)', \n",
    "         ylabel='Concentration (arbitrary units)')\n",
    "\n",
    "savefig('figs/chap18-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical solution\n",
    "\n",
    "Now let's solve the differential equation numerically using `run_ode_solver`, which is an implementation of Ralston's method.\n",
    "\n",
    "Instead of an update function, we provide a slope function that evaluates the right-hand side of the differential equations.\n",
    "\n",
    "We don't have to do the update part; the solver does it for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes derivatives of the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: derivatives of G and X\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    k1, k2, k3 = system.k1, system.k2, system.k3 \n",
    "    I, Ib, Gb = system.I, system.Ib, system.Gb\n",
    "    \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    return dGdt, dXdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run the ODE solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2, details = run_ode_solver(system, slope_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`details` is a `ModSimSeries` object with information about how the solver worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`results` is a `TimeFrame` with one row for each time step and one column for each state variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the results from `run_simulation` and `run_ode_solver`, we can see that they are not very different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(results.G, 'C0', label='run_simulation')\n",
    "plot(results2.G, 'C2--', label='run_ode_solver')\n",
    "\n",
    "decorate(xlabel='Time (min)', ylabel='Concentration (mg/dL)')\n",
    "\n",
    "savefig('figs/chap18-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differences in `G` are less than 2%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = results.G - results2.G\n",
    "percent_diff = diff / results2.G * 100\n",
    "percent_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "max(abs(percent_diff))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "**Exercise:**  Our solution to the differential equations is only approximate because we used a finite step size, `dt=2` minutes.\n",
    "\n",
    "If we make the step size smaller, we expect the solution to be more accurate.  Run the simulation with `dt=1` and compare the results.  What is the largest relative error between the two solutions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood\n",
    "\n",
    "Here's the source code for `run_ode_solver` if you'd like to know how it works.\n",
    "\n",
    "Notice that `run_ode_solver` is another name for `run_ralston`, which implements [Ralston's method](https://en.wikipedia.org/wiki/List_of_Runge–Kutta_methods)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_code(run_ode_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Related reading:** You might be interested in this article about [people making a DIY artificial pancreas](https://www.bloomberg.com/news/features/2018-08-08/the-250-biohack-that-s-revolutionizing-life-with-diabetes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
